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Abstract 

We apply the Lindstedt method to the one dimensional Fermi-Pasta-Ulam 
P lattice to find fully general solutions to the complete set of equations of 
motion. The pertubative scheme employed uses e as the expansion parameter, 
where e is the coefficient of the quartic coupling between nearest neighbors. 
We compare our non-secular perturbative solutions to numerical solutions and 
find striking agreement. 

I. NORMAL MODE REPRESENTATION 

No known solution for a coupled set of Duffing oscillators exists, with the exception 
of specific cases. The Duffing lattice is perhaps more well known as the FPU /3-system, 
one of the systems studied by Fermi, Pasta, and Ulam in the 1950s [1] to numerically 
investigate the approach to thermal equilibrium in nonlinear systems. FPU found that 
the expected equipartition of energy among the linear normal modes did not occur, and 
that periodic, non-ergodic motion persisted even in the presence of the nonlinear mixing 
terms. Subsequent investigations, both numerical and analytic, have shown that different 
parameter regions specified by the total energy and the strength of the nonlinear terms 
as specified by the parameter e lead to markedly different types of behaviour [2] , with the 
expected ergodicity emerging with increasing energy and e. But the observation of the 
periodic motion in the orginal FPU /9-system computer experiments lends credence to the 
possibility of finding acceptable pertubative solutions in certain parameter regions. That 
no exact solution exsists for a large - body lattice modulo some exceptional cases [3] 
motivates us to develop pertubative solutions to systems of Duffing oscillators. 

Several pertubation schemes as applied to the Duffing oscillator have been investigated 
[4], but the methods are applicable only in the weakly nonlinear regime and also require 
the time t to be small in order to be rigorously trustworthy. We therefore apply a differ- 
ent perturbative procedure called the Lindstedt series. We note that similar perturbation 
schemes that avoid secular terms have been proposed [5-11] in the past. In particular, we 
note that the method presented here is a more general variation of the so-called shifted 
frequency perturbation scheme, as described in [11]. The Lindstedt method is a systematic 
procedure to compute formal power series expansions of quasi-periodic solutions [12]. We 
consider an anharmonic bath of degrees of freedom in one dimension with fixed endoints 
{qo = Qn+i = 0), (we consider these boundary conditions as opposed to periodic boundary 
conditions so as to avoid zero modes, which are problematic when using the lattice as a bath 
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coupled to a test particle in the standard derivation of the generalized Langevin equation) 
as described by the Hamiltonian: 

-1 N 1 N+1 N+1 

H-^Ep' + ^E i1^ - 1^-^f + J E - ^-l)' • (1) 
i=l i=l 1=1 

The normal coordinates Qi are related to the site coordinates Qj through the transfor- 
mations 

N 

q^ = Y.^^3Q3 (2) 



and 



N 

Pi-E^ijPj (3) 



where 



N+1 \N+1 

By means of this transformation, the harmonic part of the Hamiltonian is decomposed into 
a sum of independent normal modes: 



N 1 

^0 = ^2 (^^'+^^^9' 



k=l 

where 

We have the following compact expresssion for H expressed in terms of the linear normal 
modes [13,5,14]: 



1 ^ 



N 

2 K= 

e 



k=l 

N N N N 

t 

+ 



o / AT , 1 \ E E E E ^k^l^m^nCklmnQkQlQmQm 
+ IJ fe=U=lm=ln=l 

where the coupling coefficients Ckimn are given by 

Cklmn — ^fe+i+m+n + ^fe-i+m+n + ^k+l-m+n 
+ ^fe+/+m-n + ^k-l-m+n + ^fe+/-m-n 
~l~ ^fe— i+m— n ~l~ ^k—l—m—m 

with defined by 



1 for r = 0, 

-1 for r = ±2 (AT + 1), . (7) 
otherwise 
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II. LINDSTEDT METHOD 



We seek solutions to the equations of motion for Qj. (t) , 

^2 ^iQk = ~ 2(N + 1) ^ '^k^l^m^nCklmnQlQmQn (8) 
V / l,m,n 

where e is a small parameter (e < 1). A hallmark of nonlinear, periodic behaviour is the 
amplitude dependence of the frequency. To capture this relationship, for each normal mode 
Qk we define a new stretched temporal variable = (3kt where to first order in e we define 
= 1 + epk and the quantity pk is to be determined. The Lindstedt method can in principle 
be expanded to all powers in e [12], although the convergence properties for arbitrary N and 
e are not known. In this work, we will only retain expressions up to first order in e in the 
quantities {Pk}- 

The Lindstedt method then makes the assumption that we can express Qki^k) in the 
form 

Qk ~ Qk,o + eQfc,i- (9) 
The zeroth order equations are (e = 0): 



with solutions 



+c^^Qm = o> (10) 



Qkfi{rk) = (5fe(0) cos(a;ferfe) + sin{uJkTk). (11) 

^k 



The order e equations are given by 

, 2 +^kQk,l = '^Pk^kQkfii'Tk) 



2{N + l,, 

or, in terms of the original time variable t: 



'^l^m^nCklmnQl,o{Tl)Qm,o{Tm)Qn,o{Tn) ■ (12) 



d Qk,l . o2 2/0 o /32 2 



+ (5iujiQk,i = 2pk/3iujiQk,o{M 



d'^LOk ^ 

^l^m^nCklmnQl,o{Pl^lt)Qm,o{Pm'-^mt)Qn,o{Pn'-^nt) ■ (13) 



2(N+1), ^ . 
V / l,m,n=l 

In order to apply the Lindstedt method, we need to identify all of the potential resonant 
driving terms [12]. Terms with the arguments {Picui + PmUJm + Pn^njt will never provide 
resonance terms. Using the symmetry properties of the coupling coefficients Ckimm we 
arrive at the following resonant forcing terms in the equations of motion for the Qk,i'- 
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Qk{0) cos{l3kUJkt) + sm{l3kUJkt) 



'^kkkk ^ 



8(A^ + 1) 



cos{uJkl3kt) 



+ 



sin(a;fe/3fct) 



N 



2Qk{0)Ql{0) + 2Qk{0) 



2 H 



+ 



2Q^ (0)M?) +2^^^°) ^^'"(^^^ 



+non-resonant terms. 



sm{PkUJkt) \ ■ 



(14) 



Setting to zero the coefficients mutiplying sin(/3fcC<jfct) and cos{j3kUJkt) yields the following 
system of coupled nonlinear algebraic equations for the parameters {pk}'- 



Pk = 



+ 



16(A^+ 1) 
3 ^ 



Ckkkk{ Qfc(O) + 



Qk{Q) 



{l + epk)ujk, 



XI ^IfikkmmS <3m(0) + I 77" 



gm(o) 



+ (■Pm)(^r, 



(15) 



In principle, one may solve this system of equations for the set {pk}- Having obtained 
{pk}-i one then substitutes the expressions given by Eqn. (11) into the right hand side of 
Eqn. (13) and uses the harmonic oscillator Green's function to solve for the set {Qk,i}- In 
the present work, for the sake of simplicity, we restrict our attention to corrections to the 
harmonic frequencies {ujk} that are first order in e; i.e., we consider /3k = I + epk where pk is 
independent of e. If one wishes to compute ffist order corrections to the {ook}, then e may 
be set to zero in the right hand side of Eqn. (15), and one obtains: 



Soul 



Pk = 



+ 



16(A^+ 1) 



Ckkkk { Qfc(O) + 



'Qkioy 
, ^k . 



N 



8(^+1)^. 



^m^kkmm { Qm(0) + 



and (3k = l + epk- 

Therefore the zeroth order normal mode Qk,Q{t) is modified as follows: 
Qk,o{t) = Qk{0) cos{uJkPkt) + sm{ujk(3kt) 

^kPk 



(16) 



(17) 
(18) 



where 
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2^ 



With this expression for the parameters where is independent of e, we have not 
removed resonant driving terms from Eqn. (13), but rather from the evolution equation for 
Qk,i{t) obtained by setting all of the set {(3^} equal to unity in Eqn. (13). With the {pk} 
given by Eqn. (16), we have removed resonant driving terms from the equations: 

^k ^ 



Therefore, consistency requires that we solve the Eqns. (20) to determine the Qk,i{t)- 

Having solved explicitly for the set Qk,o{t) in terms of initial conditions, we may now 
solve for the set Qk,i{t)- The equation of motion for Qk,i{t) then becomes 

^^2' ^^^^'^ ^ ~ 2(A/" + 1) ^ '^l'^m'^nCklmnQl,o{'^lt)Qmfl{'^m't)Qnfl{'^nt), (21) 
^ ' {l,m,n=l) 

where the brackets {l,m,n = 1) denote a restricted sum. All of the amplitude dependence 
of the frequencies is carried by the {Qk,o} when we only work to first order in e in the {/5fc}- 

The requirement that no secular terms appear on the right hand side, the same require- 
ment that fixed the parameters pk and the set j3k, has also guaranteed that there are no 
resonant driving terms in the above equation for Qk,i(t). This requirement of no resonant 
terms has thus transformed the orginal sum into a restricted sum. In the restricted sum, 
S{I,m,n=i>' discard all resonant driving terms. A fixed set of indices l,m,n will produce 
eight different harmonic source terms. The meaning of the restricted sum is that we only 
keep the harmonic driving terms where the arguments do not equal ±0;^. The resonant 
terms have already been accounted for in our definition of pk- 

The full solution Qk{t) is parametrized by two arbitrary constants determined from initial 
conditions, which implies that the solution for Qk,i{t) will contain no arbitrary constants. 
Therefore we are only interested in the particular solution for Qk,i{t)- 

Rewriting the equation of motion as 



+ ^lQk,i{t) 



N 



COk ^ ^1 ^...^ ...^ ... , Ql{0)Qm{0)Qn{0) 



S(N + 1'^ ^ ^l^m^nCklmn 



Ql{0)Qm{0)Qn{0) + 



{l,m,n=l) 

^ Qm(0)4(0)Qn(0) Qn(0)Qm(0)4(0) 



COs{{u}i+U}m - U}n)t) 
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+ 



Ql{0)QmiO)QniO) 



Ql{0)Qm{0)Qn{0) 



Qm{0)Ql{0)Qn{0) Qn{0)Qm{0)Ql{0) 



QliO)QmiO)QniO) - 



<5/(0)<5m(0)Qn(0) 



cos((a;/ + a;^ + a;„)t) 



^ Qm{0)Ql{0)Qn{0) ^ Qn{0)Qm{0)Ql{0) 
<3/(0)Qm(0)gn(0) 



<5K0)Qm(0)Qn(0) + 



<9m(o)gKo)gn(o) , g„(o)g^(o)4(o) 



+ 



4(0)Qm(0)gn(0) ^ Qm(0)Q«(0)4(0) 

Q/(o)g„(o)g^(o) gKo)Qm(o)Q„(o) 



4(o)Qm(o)Qn(o) ^ gm(o)gn(o)gKo) 

UJiUJrnUJn UJl 
Ql{0)Qn{0)Qm{0) Ql{0)QmiO)QniO) 



OJr. 



UJri 



g/(o)gm(o)g„(o) , g„.(o)g,(o)g;(o) 



+ 



gKo)gn(o)g^(o) , QmQmimni^) 



+ 



gKo)g„.(o)g„(o) , g„.(o)g,(o)gKo) 



+ 



^ Ql{0)Qn{0)Qm{0) ^ Ql{0)Qm{0)Qn{0) 



and using the results 



* cos(q;s) sin{ujk{t — s)) cos{uj}.t) — cos{at) 







{a + ujk){a - ujk) ' 



and 



* sin(Q;s) sn\{u)k{t — s)) sm[ujkt)a — sm.{at)ujk 



we see that the solution for Qk,iit) is given by the following expression: 



(22) 



(23) 



(24) 
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^ {l,m,n=l) y - 



Q/(0)<5m(0)gn(0) 



+ 



Qm(0)OKO)Qn(0) <?n(0)Qm(0)QKO) " 

<5Ko)(5m(o)gn(o) 



Q/(o)g„(o)g„(o) 

Qm{0)Ql{0)Qn{0) Qn{0)Qm{0)Ql{0) 



+ QmQm{0)Qn{0) - 



Ql{0)Qm{0)Qn{0) 



^ Qm(0)QKO)Qn(0) ^ Qn(0)gm(0)QKO) " 

Qi(o)<9m(o)g„(o) 



+ 



QKo)Qm(o)g„(o) + 



Qm{o)Qi{o)Qn{o) , Qn(o)g„.(o)4(o) 



+ 



cos(a;fct) - cos((a;; + o;^ + ujn)t) 



{Ui -Um-UJn+ U;k){u;i -Um-UJn- UJk) _ 



COs{u!kt) — COs((a'/ — UJm + '^n)^) 

(a;; - o;^ + + Uk) {coi -Um + ujn- ^k) _ 



+ 



4(o)gm(o)g„(o) , g„.(o)g„(o)g/(o) 



+ 



gKo)gn(o)g„^(o) gKo)g™(o)g„(o) 



+ 



OJr. 



UJk{uJi+ LOrn- i^n + ^k){^l+ ^m- ^n- ^k) 



QmQmimni^) , gm(o)g„(o)gKo) 



+ 



QmQnimmi^) gKO)gm(0)gn(0) 



+ 



0J„ 



^k{^l ~ ~ + ^k) {^l ~ ~ ^k) 



QmQnmw) , gm(o)gn(o)gKo) 



LOiUJmUn 



+ 



+ 



+ 



Sm{uJkt){uJl - LOm + LOn) - sin((cj; - CJ^ + ^n)t)uJk 
UJk{uJi - Um + Un + UJk){uJi - UJm + - '^k) 



gKo)g™(o)gn(o) , g,n(o)gn(o)gKo) 



+ 



UJlUJmUn 



^ QmQn{0)Qm{0) ^ Ql{0)Qm{0)Qn{0) 



UJr, 



uJki^i + o^m + ti^n + ^k){^i + u;^ + a;„ - uJk) 



, (25) 
(26) 



The above expression is the principal result of this work. In the following sections, we apply 
this formalism for specific N and specific initial conditions. 
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III. EXAMPLE FOR N = 2 



The major advantage of the pertubative solutions is their generahty; the full set of 
pertubative solutions to the FPU-/? chain presented above with N degrees of freedom has 
2A'" independent arbritrary constants specified by the initial conditions, as is necessary to 
represent a fully general solution. Solitons, discrete breathers, and other exotic nonlinear 
excitations that have been observed in laboratory FPU-/? -like systems may then presum- 
bably be constructed from the above general pertubative solution with appropriately chosen 
initial conditions, given that it does in fact represent the fully general solution in the small 
e limit. In this section we present some general solutions with initial conditions chosen so as 
to describe accurately these nonlinear excitations. To illustrate the method, we construct 
explicit solutions for the simplest non-trivial system, the N — 2 system. The equations of 
motion are: 



and 



Qi + Qi^-leQl-^eQiQl (27) 



Q^ + 3Q2^-leQl-^eQlQ2. (28) 



One easily finds that 

A = 1 + (g?(o) + pm) + 1^ (qUo) + ip|(o)) (29) 

P2^1 + ^e (quo) + \pm) + le (QUO) + (30) 



16 

As an example of random initial conditions, we consider the case (Qi(0) = 1/10, Pi{0) = 
1/10, (52(0) = 1 and -P2(0) = 0). We let e = 1/10. These parameter choices result in values 
for pi and p2 of 303/808 and 113/200, respectively (with dimensions of [e~^]). With these 
values for the initial conditions and e, the explicit solutions are: 

^ / X / /8303 , , , \ 400 /8303 , , , \ , , , 

QAt) = 1/10 cos sin l/6 7r U + sin sin I/GttU sin 1/6 tt ^ 

^ ' V4000 ^ ^ W 8303 V4000 ^ ' ' ^ ' " 

1 cos (i) -cos((2 V3- Ijt)) 1 

+ ^^^^^^ cos(t) 



2400 _2 + 2 V3 320000 

cos (3 1) 



320000 

1 (cos(t) -cos((2V3 + l)t)) V3 1 \/3(cos(t) -cos((l-2V3)t 



1600 2v^ + 2 4800 2 - 2^3 

1 (sin(t) (2 V3- 1) - sin ((2 V3- l)i)) V3 



2400 _2 + 2 V3 

^ V3(sinW(l-2v^)-s.n((l-2V3)0) _ ^ ^ 

4800 2 - 2 \/3 320000 ^ ' 320000 ^ ' 



1 (sin(i) (2V3 + 1) -sin((2V3 + l)i)) V3 
1600 2v^ + 2 



(31) 
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QJt) = COS f sin (1/3 7i)t] -— cos (Vst) + — cos (s Vst) 
' VIOOO ^ ^ V 640 ^ 640 \ ^ J 



1 V3(sin(V3i) (2- Vs) -sin ((2-^3)^)^3) 



12000 

1 V3l 


[^sin 1 




2 - 

(v^-2) 


2^3 

1 - sill ^ 








24000 

1 v^(^ 
8000 


-2 + 2v^ 

dn (VSi) + 2) - sin ((v^ + 2) t) Vs) 
2a/3 + 2 



(32) 



In Fig. (1), we plot Qiii) and Q2{t) and in Fig. (2), we plot the numeric solutions for 
the same system. 
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FIGURES 



Lindstedt 'Solutions 




FIG. 1. Qi{t),Q2{t) for the initial conditions given in the text. Qi{t) is represented by the 
lower-ampHtude curve. 



NumericarSolutions 




10 



FIG. 2. Numerical solutions for Qi{t) and Q2{t) for the initial conditions given in the text. 
Qi{t) is represented by the lower-amplitude curve. 

In Figs. (3,4, 5), we plot both the numerical solution and the analytical approximation 
for Qi for three different time intervals of twenty units, and In Figs. (6, 7, 8), we plot both 
the numerical solution and the analytical approximation for Q2 for the same time intervals. 

Q[1] 




FIG. 3. Numerical solutions for Qi{t) and the analytical approximation for Qi{t) for the initial 
conditions given in the text. 
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Q[1] 




FIG. 4. Numerical solutions for Qi{t) and the analytical approximation for Qi{t) for the initial 
conditions given in the text. 

Q[1] 




FIG. 5. Numerical solutions for Qi{t) and the analytical approximation for Qi{t) for the initial 
conditions given in the text. 
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Q[2] 




FIG. 6. Numerical solutions for (52 (i) and the analytical approximation for Q2{t) for the initial 
conditions given in the text. 



Q[2] 



0.5 



2^1 




20E 



21( 



215 



220 



-0.5 



FIG. 7. Numerical solutions for Q2{i) and the analytical approximation for Q2{t) for the initial 
conditions given in the text. 
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Q[2] 




FIG. 8. Numerical solutions for (52 (i) and the analytical approximation for Q2{t) for the initial 
conditions given in the text. 

In Fig. (9) we plot 

Q2,o(r2 = m + eQ2,i{t) - Q2,num(t) (33) 

and Q2,num(^)5 , where Q2,num(i^) is the numerical solution for Q2(t) of the system of 
eqns. (27,28). 

In Fig. (10), we plot the quantities Q2,mim(t) (black curve), our perturbative approxi- 
mation given by Eqn. (32) (grey curve), and the their difference, Eqn. (33). From Fig. (10), 
we note that the principal source of the growing disagreement between our approximate 
solution and the numerical solution arises from the difference in phases between the two 
functions. The amplitudes appear to be in strong agreement. Furthermore, from Fig. (9), 
it appears that the phase difference is some linear function of time. We thus speculate 
that the full numerical solution has not only an amplitude-dependent frequency, but also 
a time-dependent frequency. If we allow for time-dependent frequency shifts pk, we may 
achieve greater agreement for much longer times. We leave this open question for a further 
investigation. 
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FIG. 9. Plots of Q2,num(i) (thin curve) and Eqn. (33) (thick curve) . 




FIG. 10. Plots of {Q2. num(^) (thick curve), our perturbative approximation given by Eqn. (32) 
(thin curve), and the their difference, Eqn. (33). 
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IV. CONCLUSION 



The FPU-/3 lattice has been a subject of study in relation to fundamental issues in non- 
linear dynamical systems since its introduction. Because a single Duffing oscillator admits 
of an exact solution, one may be tempted to entertain the notion that a lattice of Duffing os- 
cillators with nearest neighbor couphng may also allow fully general analytic solutions. For 
special initial conditions, the full lattice dynamics is described by a single quartic oscillator 
and thus admits exact analytic solutions [15,16], but these particular solutions correspond to 
special initial conditions. While for a small number N of degrees of freedom, some of these 
special solutions are stable over an appreciable range of e, they are actually unstable for 
arbitrarily small e in the limit of approaching infinity. For more general initial conditions, 
we have applied the Lindstcdt method to lowest order in the nonlinearity parameter e to 
obtain analytical, quasiperiodic approximations. These quasiperiodic solutions are predicted 
to exist on the basis of the KAM (Kolmogorov, Arnold and Moser) theorem and is the ac- 
cepted explanation for the nonchaotic behavior and recurrences that puzzled Fermi, Pasta 
and Ulam. As noted above, the formal perturbative scheme presented above has unknown 
convergence properties. For small e and early times, the solutions of the formalism and nu- 
merically obtained solutions are numerically indistinguishable. For the example in the text, 
the pertubative and numerical solutions begin to noticeably differ in phase for dimensionless 
time on the order of forty. We conclude that the formal perturbation expansions developed 
here, truncated at first order in e, give useful analytic approximations to quasi-periodic so- 
lutions of the FPU-/9 lattice for moderately long times. Finally, we remark that the general 
application of the above perturbative scheme is independent of the spatial dimension of the 
system. 
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